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In this paper, the implementation of two wall function models in the Reynolds averaged 
Navier-Stokes (RANS) computational fluid dynamics (CFD) code FUN3D is described. 
FUN3D is a node centered method for solving the three-dimensional Navier-Stokes equa- 
tions on unstructured computational grids. The first wall function model, based on the 
work of Knopp et al., is used in conjunction with the one-equation turbulence model of 
Spalart-Allmaras. The second wall function model, also based on the work of Knopp, is 
used in conjunction with the two-equation k-ij turbulence model of Menter. The wall func- 
tion models compute the wall momentum and energy flux, which are used to weakly enforce 
the wall velocity and pressure flux boundary conditions in the mean flow momentum and 
energy equations. These wall conditions are implemented in an implicit form where the 
contribution of the wall function model to the Jacobian are also included. The boundary 
conditions of the turbulence transport equations are enforced explicitly (strongly) on all 
solid boundaries. The use of the wall function models is demonstrated on four test cases: a 
flat plate boundary layer, a subsonic diffuser, a 2D airfoil, and a 3D semi-span wing. Where 
possible, different near-wall viscous spacing tactics are examined. Iterative residual con- 
vergence was obtained in most cases. Solution results are compared with theoretical and 
experimental data for several variations of grid spacing. In general, very good comparisons 
with data were achieved. 
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I. Introduction 


The increasing size of numerical simulations of steady and unsteady aerodynamic flows demands faster 
and more efficient computational methods to better use available resources. Robust and compact numerical 
schemes can reduce computational time and memory requirements. A reduction in computing requirements 
can also be achieved through solving a smaller problem by using models that make simplifying physical 
assumptions that reduce the grid resolution requirements in certain regions of the flow. 

Historically, the log-law of the wall for turbulent flows has been used to construct a wall function model, 
which is then used to significantly reduce the grid resolution requirements in the immediate vicinity of walls. 
This reduction in the required wall grid resolution is achieved because the first grid point off the wall can 
be located in the log-law region of the boundary layer, 25 < y+ < 250, instead of the viscous sub-layer, 
y’*' <= 1, where y"*" is a scaled wall-normal coordinate or wall units. The wall function is used to provide the 
wall momentum flux and wall energy flux to the Navier-Stokes equations; and can be used in conjunction 
with any number of turbulence modeling strategies from linear eddy viscosity models, hybrid RANS/LES^ 
or DES.2 

However, the practical application of log-law wall function models to complex flows has some challenges, 
the most notable of which is the difficulty of knowing where to place the first grid point off the wall when 
generating the grid. In complex flows and/or geometries, the state of the boundary layer is not known a priori 
making it difficult to generate a grid such that a log-law wall function treatment is well posed everywhere 
in the computational domain. A more practical approach would be to adapt the grid such that the first 
point off the wall is always located at a y+ that makes the use of the wall function well posed. However, 
another shortcoming of log-law wall function models is that they are not valid when the wall shear stress 
becomes zero, for example points of flow separation or re-attachment. Therefore a wall function model that 
adapts to the physics of the flow local to the wall is desired. Esch et al.,^ Menter et al.,"^ Kalitzin et al.® and 
Knopp et al.® have developed and demonstrated wall function models that blend viscous sub-layer function 
models and log-law function models together such that the resulting hybrid wall function model is more 
universal and thus, well posed for practical values of y+, e.g., 0 < y"*" < 250. Moreover, Knopp et al.® have 
demonstrated the use of this approach in conjunction with a near wall grid adaptation scheme. 

The objective of this paper is to describe the implementation of the hybrid wall function models of Knopp 
et al.® in the node-based unstructured-grid CFD code Fun 3D and compare their accuracy relative to theory 
and relative to solutions without a wall function model for several flows of interest. The first section will 
discuss the turbulence and wall function models, followed by the implementation of the wall function in 
Fun 3D, test case descriptions and a discussion of the results. 

II. Turbulence and wall models 

The turbulence models of Spalart-Allmaras^’ ® and Menter et ah,® are used in conjunction with the wall 
function approach published by Knopp et al.® Both of these turbulence models are well established with 
the aerospace community and have been used to simulate flow around a wide variety of geometries and flow 
situations. 


A. Spalart-Allmaras One-Equation Model 

The original one-equation model of Spalart-Allmaras is given by Eq. 
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where the turbulent eddy viscosity is 
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This study also used a modified form of the Spalart-Allmaras model. The modified model is the same as 
the original version except for when the turbulence variable v is negative. In the situation of non-positive 
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values, Eq. I^is solved instead. 

— P + V(t/i>) =Cf,i(l-Ct3)|W|i> + Cu,i f + - [V • ((i^+ C’/„) Vi>) + Cfe2Vi> • V;>] (3) 
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The resulting turbulent eddy viscosity (/it) is set to zero when i> is negative (i> itself becomes a passive 
scalar). Other model parameters are defined as follows: 
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The constants are: 
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Cyj2 — 0.3, c,„3 — 2, Cyi — 7.1, Ct3 — 1.2, Ct4 — 0.5 

Cbl , 1 + Cb2 

^wl — y “t“ 


( 8 ) 

(9) 

(10) 


B. Menter-SST Two-equation model 


The two-equation model used in this paper is the model proposed by Menter.® Eqs. 11 and 12 are the 
transport equations written in conserved variable form. 
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The production of turbulent kinetic energy is stress-based 

Pk = min (Pk, 10 (3* pujk) 
with the turbulent eddy viscosity computed from: 
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Each of the constants, ct^, 7 and /3, are a blend of an inner (1) and outer (2) constant, blended via: 
(j> = Fi(j>i + {1 — Fi)(j )2 where t/i represents constant 1 and (p 2 represents constant 2. Additional functions 
are given by: 


El = tanh (r^) , Pi = min 


max 


Vk 500 j/\ 4pCToj2k 

/3W ] ’ CDkV? 


( 16 ) 


3 of[27] 


American Institute of Aeronautics and Astronautics 


CDktj = max 


(17) 
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The turbulent variable boundary conditions are: 


The constants are: 
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Further discussions of these models can be found 


on the turbulence model resource websited^’^^ 
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C. Wall function model 


The wall function model proposed by Knopp et al.® involves two general steps to compute the friction velocity, 
Ur|wt- In the first step, transcendental solutions for the friction velocity using the canonical equation of the 
log-law, Eq. 19 and the law-of-the-wall equation by Reichardt,^^ Eq. 20 are blended using a hyperbolic 
tangent function, Eq. |21[ 

Equations [T^ and |20| are shown in Fig.[T] as the dashed and dashed-dot curves respectively. As reference, 
the horizontal brace shows possible wall spacings for use with the wall function model, ranging from less 
than one to approximately 500. We see that Reichardt’s equation already matches the viscous sub-layer 
quite closely, as seen by the overlap of the solid and the dashed-dot curves. 


u^ = y^ 

u" = ln(y")/K + 5.1 
Reichardt 



Figure 1: Canonical boundary layer profiles. 


I Log 


I Reichardt 


-In (y(^) -f 5.1. 

ln(l-f 0.4y+) 

— ^ -^ + 7.8 


1 _ ^ 


11 



(19) 

( 20 ) 



American Institute of Aeronautics and Astronautics 


The resulting blended model, Reichardt,m, shown as the dashed-dot curve in Fig. Ha), overlays both the 
sub-layer and the log-layer regions. For reference, the hyperbolic tangent blending function, (pbi, is plotted 
along side the blended wall function, in the same figure. 


u" = ln(y")/K + 5.1 

Reichardt 

Reichardt,m 



Reichardt, m 
Spalding, 3 terms 
Blended Spalding, SST 



Cb 

Cf) 

8- 


(a) Blended Reichardt step. 


(b) Blended Spalding step. 


Figure 2: Schematic of blended profiles. 
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In the following step, the friction velocity from Eq. m and the friction velocity from the transcendental 
solution of the Spalding friction velocity, Eq. 22 or Eq. 24 are blended using a hyperbolic tangent functions 
Eq. or Eq. The wall model, as applied to the SA model, uses a 5-term expansion, while as applied 
to the SST model, a 3-term expansion is used. Also, slightly different functional forms of the blending are 
used depending on the turbulence model. 

For SA: 


rr-T I Spalding, 5 ^ Yl — 

MtIsA = (1 - < 
For SST: 


+ 


^ — k5.2 / „KU 


_ 1 _ KU+ - 


(ku)^)^ (ku)^)^ (ku)^)"^ (r«Ui) 


6 


24 


I Spalding, 3 


y+ = U+ -b - 1 - KU+ - 


+ (^u+)^ (rcu+)3) 


6 


WtIsST = (1 — <PSSt) Wt I Spalding, 3 + y^SST^^r |Roicliardt,m, V^SST = tank (arg^) , arg = 
The boundary condition for uj, Eqn. 34 in Knopp,® is computed as 

y + 
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D. Compressibility 

The effects of compressibility are modeled via a correction to the wall tangential velocity, ui, by using van 
Driest II transformation,!® 
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where is used in place of ui in any transcendental equation involving the solution of the friction velocity 
at the wall. 


III. Computational Method and Implementation 


A. Fun 3D Code 

FunSD is an unstructured three-dimensional, implicit, Navier-Stokes code. Roe’s flux difference splitting!® 
is used for the calculation of the inviscid terms. Other available flux construction methods include HLLC,!!" 
AUFS,!® and LDFSS.!® The default method for calculation of the Jacobian is the ffirx function of van 
Leer,^° but the method by Roe and the HLLC, AUFS and LDFSS methods are also available. The use of 
flux limiters are grid and flow dependent. Flux limiting options include MinMod^! and methods by Barth 
and Jespersen^^ and Venkatakrishnan.^® Other details regarding FunSD can be found in Anderson and 
Bonhaus^! and Anderson et al.,^® as well as in the extensive bibliography that is accessible at the FunSD 
Web site.^® 

B. Boundary Element Discretization 

Discretization of the computational volume consists of any combination of 4 basic polyhedra: tetrahedral, 
hexahedral, prismatic, or pyramidal elements. The faces of the volume elements are either triangles or 
quadrilaterals. The boundaries, therefore, consist only of triangular or quadrilateral faces or combinations 
of both types of face elements. Each of these faces are divided in to three or four sub-elements, or duals, 
depending upon being triangular or quadrilateral, respectively. 


z z 



Figure S: Schematic of a wall boundary element configuration. 
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external flow 



Figure 4: Schematic of domain decomposition (Figure adapted by author from® with permission). 


Each node on a boundary face is associated with a number of quadrilaterally shaped sub-elements, or 
duals. For example, four boundary duals, Di through D 4 , associated with the node Pi, are sketched in 
Figure]^ a), where Pi, P 2 , P 3 , P 5 , Pe, and P 7 are all boundary nodes. 

Figure. ib) depicts a pyramidal element with one face associated with the boundary. Several other 
neighboring tetrahedra are sketched in with lighter weight lines. The heavier black lines, defined by Pi, P 2 , 
P 4 , Pe, and P 7 , are the edges of the example pyramid. Pi, P 2 , Pe, and P 7 are the boundary nodes of the 
example pyramid. (A shorthand description for this collection of nodes will be [1,2, 6 , 7].) 

Four boundary elements have an edge associated with Pi, one quadrilateral (quad) [1,2, 6 , 7], and three 
triangles (tri) [1,5,6], [1,3,5] and [1,2,3]. The green quad, Di, is one of the previously mentioned duals on 
the boundary face of the quad and would be used as part of the the construction of the flux due to the 
wall function model. The vector, di, represents the minimum spacing value associated with the boundary 
node Pi. The minimum spacing is the average distance of all the off-boundary nodes of the element. In this 
example with the pyramid, it is the distance off the wall of the point P 4 . 

The sample pyramidal element could occupy almost any region of a wall-bounded flow, schematically 
shown in Fig. depending upon conditions and geometry. In this example, the element extends from the 
wall, Fw, to some point in the log-layer, F^. The region below F^ is modeled by the wall function and the 
outer region, fl, is resolved by the general flow equations. Wall function blending allows accommodation of 
grid spacing spanning from the viscous sub-layer to approximately 250 wall units above the solid boundary. 

Depending upon the type, boundary conditions in Fun3D are applied at either the boundary nodes, e.g.. 
Pi, or across the boundary element duals, e.g., Di. Strong boundary conditions, such as the viscous no-slip 
wall boundary, in Fun3D nomenclature the viscous_solid boundary, explicitly set the velocity at all the 
boundary nodes. Pi, to zero. The values of the turbulence model variables, v, k, and u are also explicitly 
set at the wall boundary nodes. A weak boundary condition, such as the viscous no-slip wall boundary 
and the wall function boundary, viscous_weak_wall and viscous_wall_function in Fun3D nomenclature, 
respectively, is applied implicitly at Di. The viscous stresses and pressure fluxes are determined for all the 
boundary duals and added to the right hand side of each of the respective momentum equations. Additionally, 
to enforce a no-slip velocity condition, a zero momentum flux is assumed across the boundary element. 



American Institute of Aeronautics and Astronautics 


C. Wall Function Implementation 

The wall function and the weak wall boundary condition influence the flow through the stress and pressure 
terms in the momentum equations. The pressure term is computed and applied in the same way between 
the two boundary types. With a wall function, the shear stress is modeled using a number of transcendental 
equations. The weak wall boundary condition, by comparison, computes a shear stress term directly from the 
flowfleld solution and grid element geometry. This section describes the wall function boundary condition. 


Solution process 

( 1 ) For each element attached to a wall function boundary, an averaged density, velocity, pressure and 
viscosity is computed. Just as the minimum spacing value is determined by the average distance of 
all the off-boundary nodes in the element, the primitive flow variable values are the average of all the 
off-boundary values and are also designated by the subscript “1”. 


qi 


1 


nodesoi 
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nodesoff “ 

n—1 


q= [p, 


sp] 


(29) 


( 2 ) The off-wall velocities determine the local flow direction from which a tangential velocity is computed. 

U|| = ui — u_Lri, u_L = ui • h, ui = [m, v, w] (30) 

1 

The unit direction of the flow tangential to the wall is then 

- '^11 


(31) 


( 3 ) The temperature and density are adjusted according to the magnitude of the tangential velocity. 


P = Pt^, T = T- i(7- 1) Prt |u||P, 
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and Sutherland’s law is used to determine the dynamic viscosity. 

C* = 110.3/Tref, T„f = 288. 15K. 
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The friction velocity, Ut, is computed solving Eqs. ^ 20 and 21 simultaneously via Newton’s method. 


If the distance from the boundary in wall units, yj^ = (rtrdi/f'), is less than 250, Newton’s method is 
again used, this time to solve Spalding’s relationship, Eq. along with blending equation, Eq. 
when using the SA turbulence model. Eqs. [M| and [25] are used for the SST turbulence model. 


If the distance from the boundary is greater than 250, Spalding’s relationship becomes ill-conditioned 
and difficult to solve. In that case, step (5) is not executed, and the friction velocity from step (4) is 
used. 


(6) The shear stress on the boundary is 


^wf — P^rlwf 


(34) 


(7) The shear stress, Twf, computed from the wall function, the pressure, p, and the area of the element dual 
on the wall. A, are used to compute the momentum flux contribution. Heat transfer is not accounted 
for in the current formulation. 


= rwfU||A-pnA 


( 35 ) 
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IV. Test Cases 


The four geometries used to validate the implementation of the wall function models are described in 
sub-sections through [P] Table lists relevant aspects of the simulation for each case. The flat plate and 
asymmetric diffuser flows are low speed with little, if any, compressibility effects. The asymmetric diffuser 
flow has both flow separation and re-attachment points. The flow about the RAE 2282 airfoil is attached 
for the specific flow conditions chosen for this report, and the flow about the ONERA M6 wing is a 3-D 
problem with flow separation. Comparisons are made with historic, canonical data or experimental results 
where available. 


Table 1: Test geometry setup. 


Geometry 

Grid 

Elow physics 

Flat plate 

2D 

incompressible, dp/dx « 0, affached flow 

Asymmetric subsonic diffuser 

2D 

incompressible, separafed flow 

RAE 2822 airfoil 

2D 

compressible, affached flow 

ONERA M6 wing 

3D 

compressible, separafed flow 


A. Zero pressure gradient flat plate 

In this section, the initial validation tests of the wall function model are discussed. Friction velocity, local 
skin friction and mean flow profiles are compared with results from a wall integration turbulence model for 
the flat plate and are shown in Figs. 0111 and 1^ respectively. The determination of the friction velocity was 
explored using three approaches for specifying grid spacing and grid density: 

1. Uniform refinement 

2. Variable wall spacing with fixed grid dimensions, log-layer unpacked (LLU) 

3. Viscous layer node removal from a baseline grid, log-layer packed (LLP) 


0.0003 r 


0.0002 


0.0001 - 


0 i- 


V 

Uniform 
refinement, 
y; « 0.4 ~ 3 




LLU, 0-LLP, 

y ;« 1~8 y ;« 1-^11 


Figure 5: Comparison of near-wall node spacing approaches. 


A sample of node distributions for these approaches are sketched in Fig. The generation of uniformly refined 
grids is used for assessing grid convergence order properties, see the first four lines on the left. Producing 
wall function grids through successive coarsening of a very fine baseline, though, can inadvertently reduce the 
resolution in other regions of the flow. Simply increasing the wall spacing, referred to as log-layer unpacked^^ 
(LLU) spacing, can produce grids with the same issues as successive coarsening. These grids, though, are 
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often more convenient to generate. The middle four lines in the figure are examples of LLU node spacing. 
Counting from the left, the fourth and the seventh lines appear to have quite similar distributions. Node 
removed grids, or LLP spacing, maintain resolution away from the wall, as shown with the last four lines 
on the right side of the figure, while providing a large initial wall distance. The particular LLP examples 
shown use the finest uniform refinement grid as the starting point and successively remove a number of nodes 
between the wall and the outer region. 


Table 2: Grid parameters for the flat plate. 


Level 

Dimensions 

(baseline) 

Nodes 

(baseline) 

Nodes 

removed 

Spacing 
range [m] 

range 

0 

897 X 449 

402753 


- 

4.0 xlO"® 

0.4 

1 

449 X 225 

101025 


- 

8.0 xlO"® 

0.8 

2 

225 X 113 

25425 


- 

1.6 xlO"® 

1.5 

3 

113 X 57 

6441 


- 

3.4 xlO"® 

3.2 

LLU 

257x129 

33153 


- 

1.16 xlO"® - 1.05 X 10"® 

1. - 98. 

0-LLP 

897 X 449 

402753 

0' 

-40 

4.0 xlO"® - 8.4 X 10"'‘ 

0.4 - 23. 

1-LLP 

449 X 225 

101025 

0' 

-40 

8.0 xlO"® - 8.4 X 10-4 

0.8 - 76. 


The parameters for each of the gridding approaches are listed in Table The wall spacing of the uniform 
refinement series spanned from roughly one-half to 3 wall units. Each level removes every other node in 
the streamwise and normal directions. Level 0 is the finest grid. With every 8 nodes removed. Level 3 is 
the most coarse grid. The LLU approach fixed both the streamwise and vertical number of nodes in the 
grid and only increased the initial wall spacing, and then continued expanding away the wall. As mentioned 
previously, grid density near the wall is reduced with this method, possibly compromising resolution of flow 
features in the region. The baseline grid for the LLU series was between the spacing of the level 1 and 2 
uniform refinement series. A total of 13 grids were generated for the LLU series with the parameters listed 
in Table 121 

A total of 11 grids were generated for the each of the two log-layer packed (LLP) series. Between 0 and 
40 layers of nodes were removed from two of the uniform refinement grids. This method, therefore, retained 
the same distribution of nodes in the region above the wall function space. The 0-LLP series, using the level 
0 uniform refinement grid as the baseline, spanned the wall spacing range from 0.4 to 23, depending upon 
the number of layers removed. The 1-LLP series, using the level 1 uniform refinement grid as the baseline, 
spanned the wall spacing range from 0.8 to 76, again, depending upon the number of layers removed. The 
grid parameters are listed in Table All of the grids consisted of hexahedral elements and were one-cell 
wide. 

The 1-equation turbulence model of Spalart-Allmaras^’ ® and the 2-equation turbulence model of Menter® 
were used. The simulation conditions were Moo = 0.097 and ReL = 2.7 x 10®, (L = 1. m). All solutions were 
run second-order spatially accurate and in 2-dimensional mode. Solutions were obtained on the uniformly 
refined grids using both the wall integration and wall function models. Only wall function solutions were 
obtained on the LLU and LLP grids. 

Iterative convergence was obtained for all solutions with the grids listed in Table |2] Each solution 
residual history is not shown here, but the set of continuity equation residual histories shown in Fig. are 
representative of the convergence obtains for all the solutions obtained in this section. The Spalart-Allmaras 
turbulence model converged within, 1000 or so iterations, while on the order of 10,000 iterations were required 
to obtain a similar level of convergence using the SST turbulence model. With model validation being the 
primary intent of this study, more effort was extended towards obtaining iterative convergence than finding 
an optimal solution path. 

The first test examines the convergence of the predicted friction velocity with grid size. When applied 
to the uniformly refined grids, the wall function predicted the same friction velocity as the wall integration 
models, as shown by a comparison of the circle and square symbols in Fig. regardless of turbulence model. 
Identical convergence, for this sequence of grids, indicates that solutions can be obtained on wall integration 
grids with wall functions. 
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(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 6: Solution residual, continuity equation, flat plate. Moo = 0.097, Ren = 2.7 x 10® 
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(a) SA turbulence model. 
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(b) SST turbulence model. 


Figure 7: Comparison of predicted friction velocity, x = 4.0m, Moo = 0.097, RCx = 1.1 x 10^ 
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Wall integration, y* = 0.8 
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(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 8: Comparison of predicted mean flow profiles, x = 4.0m, Moo = 0.097, Rex = 1.1 x 10^ 
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(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 9: Comparison of predicted local skin friction. Moo = 0.097, Reu = 2.7 x 10® 
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The first three grids for the LLU series, in the range 1 ^ yj*" ^ aligned closely with the uniform 
refinement grids. The larger initial wall spacing results diverged in the region of the buffer for the SA 
results, and generally above the buffer region for the SST results. Results obtained on LLU grids were 
inconsistent at best. 

Nearly constant friction velocities were predicted by solutions obtain on the series of LLP grids using 
the SA based wall function. Solutions based on the SST wall function were less constant across the viscous 
sub-layer, compared to the SA solutions. The predicted friction velocity did finally level out when the initial 
spacing was above yj*" = 8. The SST model varied as much as 2%, from the finest to the coarsest, and was 
more consistent than the results using the LLU grids; as shown by a comparison of the diamond to the delta 
and triangle symbols in Fig. 

Computed wall integration and wall function mean flow velocity profiles using the LLU grid are plotted in 
wall units and compared with empirical data in Fig.[^ The wall integration solution is from the level 1 refined 
grid. Though the calculations using the LLU grid show a relatively wide variation in the predicted friction 
velocity, the mean flow velocity profiles were very close to the wall integration profile. The calculations using 
the other uniform refinement and LLP grids are not plotted but were similar to these results. 

Computed wall integration and wall function local skin friction distributions using the LLU grid are 
plotted with empirical data in Fig.[^ Similar to the mean flow profiles plotted previously, the wall integration 
solution is from the level 1 refined grid. Downstream of the fixed transition, the LLU local skin friction closely 
matched the wall integration values. The calculations using the other uniform refinement and LLP grids are 
not plotted but were similar to these results. 

B. Asymmetric subsonic diffuser 

The second test case was the asymmetric diffuser of Buice and Eaton. The adverse pressure gradient due 
to the diffuser is sufficient to cause the flow to separate. The flow subsequently re-attaches downstream of 
the diffuser in the constant area outflow section. Experimental mean flow and shear stress profiles, as well 
as surface static pressure and local skin friction data, are compared with solutions obtained with both wall 
integration and wall function models. 

The reference Mach number was set at 0.06. The Reynolds number per reference length was 20,000, 
where the reference length was the upstream channel height, H. The diffuser section was 77 H long, preceded 
by a 110 H long channel. The wall integration grid was a single-block, single-cell wide, hexahedral type grid. 
Two wall function, LLP, grids were created by removing the first 10 and first 20 wall layers, respectively, 
from the baseline wall integration grid. Layers were removed from both the upper and lower walls. Figure [T0| 
is a sketch of the diffuser geometry, along with detail views of the start and end of the diffuser section for 
the 20 node removed grid. 

Iterative convergence was obtained for all solutions. Continuity equation residuals are shown in Fig. 
for all the grid and turbulence model variations. With model validation being the primary intent of t 
study, more effort was put towards obtaining iterative convergence than finding optimal solution paths. 

Local wall spacing distributions are shown in Fig. for each of the grids. The average spacing of the 
upper wall was similar to the lower wall, so for clarity, only the lower wall distributions are plotted in the 
figure. In the region of the diffuser, 0 < x/H < 20, the first cell spacing in wall units varied from an average 
of 1 to 20, depending upon location. The local minima in y]'’ around x/H = 4 and 30 are indicative of the 
flow separation and re-attachment points. Generally the same distribution of y/" for each grid type occurred, 
regardless of turbulence and wall function model. 

Figure shows the velocity profiles in the diffuser with and without wall functions. For purposes of 
visualization, the vertical and horizontal axes are scaled differently so that the velocity profiles are shown in 
the context of the diffuser geometry. The velocity is normalized by the bulk velocity, Ub, at the entrance to 
the diffuser, x/H = 0. Computed turbulent shear stress for each wall model is compared with experimental 
data in Fig. A similar visualization is used here where the vertical and horizontal axes are scaled 
differently so that the shear stress profiles are shown in the context of the diffuser geometry. The shear 
stress is normalized by the square of the bulk velocity, Ub. Generally, the same result was achieved across 
grid types and turbulence models. 

Computed local skin friction on the lower wall of the diffuser, for each wall model, is plotted with 
experimental data in Fig. |15| Solutions obtained with the SST model separated from the diffuser wall 
further upstream of the experimental data than solutions obtained using the SA model. It is fortuitous 
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Figure 10: Sketch of asymmetric subsonic diffuser. 
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(a) SA turbulence model. 
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Wall Integration, baseline 
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(b) SST turbulence model. 


Figure 11: Solution residual, continuity equation, asymmetric diffuser. 
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Wall Integration, baseline 
Wall function, baseline 
Wall function, 10 nodes removed, LLP 
Wall function, 20 nodes removed, LLP 


Wall Integration, baseline 
Wall function, baseline 
Wall function, 10 nodes removed, LLP 
Wall function, 20 nodes removed, LLP 




x/H 


x/H 


(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 12: Upper and lower wall spacing distribution, asymmetric subsonic diffuser, Mqo = 0.06, Re = 

20,000. 


o Data, Bulce and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 

Wall function, 20 nodes removed, LLP 



x/H + 10 U/U^ 

(a) SA turbulence model. 


o Data, Bulce and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 

Wall function, 20 nodes removed, LLP 



x/H + 10 U/U^ 

(b) SST turbulence model. 


Figure 13: Velocity profiles, asymmetric subsonic diffuser, M^o — 0.06, Re = 20,000. 
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Data, Buice and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 

Wall function, 20 nodes removed, LLP 


Data, Buice and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 

Wall function, 20 nodes removed, LLP 





x/H + 300 uV/U^ x/H + 300 uV/U^ 


(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 14: Turbulent shear stress profiles, asymmetric subsonic diffuser, Moo = 0.06, Re = 20,000. 


o Data, lower wall, Buice and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 



o Data, lower wall, Buice and Eaton 

Wall Integration, baseline 

Wall function, baseline 

Wall function, 10 nodes removed, LLP 
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x/H 


(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 15: Local skin friction, asymmetric subsonic diffuser, Moo = 0.06, Re = 20,000. 
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that the point of flow separation computed using the 20 node removed grid is closer to the data than the 
computation with 10 node removed grid. 

C. RAE 2822 airfoil 

The third test case is the RAE 2822 airfoil.^® Wall integration and wall function solutions were obtained 
on a set of structured LLU grids and a separate set of structured LLP grids, which were converted to an 
unstructured format for use in this study. The first set of solutions were obtained on a sequence of grids with 
constant node count, with wall spacings ranging from approximately 1 for wall integration calculations to 10, 
20 and 40 wall units (on average) for the wall function calculations. The second hexahedral grid came from 
the NPARC Alliance validation archive. In addition to the baseline grid, 5 layers of nodes were removed 
from the grid layers adjacent to the airfoil surface. Solutions using the “node removed” grid, denoted again 
with LLP, are shown along with the initial set of LLU grids. 

Distributions of grid, wall spacing (yj*"), surface pressure coefflcients, solution residual and skin friction 
are shown in Figs . [16| through A hyperbolic grid generator was used to create the first set of C-type grids. 
All of the grids had 161 grid points normal to the airfoil surface, 513 points on the airfoil and 64 points in the 
wake region (upper and lower), for a total of 641 points in the tangential direction. For the yj*" ~ = 1 grid, 
the distance of the first point from the surface is 0.000005 chords. The step size in the tangential direction at 
the leading edge is 0.0008 and at the trailing edge is 0.0012 chords. The far-fleld and downstream boundaries 
were approximately 20 chords away. The geometric progression ratio near the surface was approximately 
1.1. The wall function grids were LLU type, so they were created by merely increasing the spacing at the 
wall, as listed in Table All other grid generation parameters were unchanged. As stated earlier, the 
second hexahedral grid came from the NPARC Alliance validation archive with the parameters also listed 
in Table m 


Table 3: Grid parameters, RAE 2822 airfoil. 


Grid 

Normal spacing [chord] 

yf 

Dimensions 

Nodes 

Wall integration 

5.0 xlO"® 

1 . 

641 X 161 

103201 

Wall function (LLU) 

5.0 xlO"'^ 

10. 

641 X 161 

103201 

Wall function (LLU) 

1.0 xlO"'^ 

20. 

641 X 161 

103201 

Wall function (LLU) 

2.0 xlO"'^ 

40. 

641 X 161 

103201 

Wall integration (NPARG) 

1.2 xlO"'^ 

2.5 

369 X 65 

47904 

Wall function (LLP) 

1.3 xlO"'^ 

25. 

369 X 60 

44214 


A compar ison of the distribution of computational nodes normal to the airfoil around x/c 
Both grids have roughly the same wall spacing of yj*" 


shown in Fig. 


16 


0.25, is 

20 to 25. Whether or not this kind 
of difference in grid distribution makes a difference is likely quite problem dependent. Iterative convergence 
was achieved for all grid variations for both the SA and SST turbulence model/wall function combinations 
run, see Fig. pT} 

Wall spacing, (yj^), distributions for all of the grids are plotted in Fig. 18 No effort is made to differentiate 
upper surface values from from lower surface values, as the intent of the plot is merely a general evaluation 
of the wall spacing trend. Exact averages were not calculated, but in general, the wall spacing annotation 
for each grid is a reasonable representation of the actual spacing of the calculation. 

Surface static pressure distributions for all the grid variations are plotted in Fig. |19| No significant 
differences in the solutions are seen. Upper surface local skin friction distributions for all the grid variations 
are plotted in Fig. The computed skin friction distributions from the six different grids were very close 
in level and variation. It is difficult to draw a significant conclusion from these plots alone. The integrated 
forces and moments for each of the grids are listed in Tables and The solution obtained using the wall 
function (LLP) grid had the largest deviation, using the wall integration yj*" = 1 as the basis for comparison, 
being off by 3% in lift, 2% in drag and 5% in pitching moment. The lift, drag and pitching moment, obtained 
with the LLP grid, was below the results from the yj*" = 20 and yj*" = 40 grids for lift, drag and moment. 
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0.058 

N 
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0.247 0.248 0.249 0.25 0.251 


X 


(a) yj*" ~ 20, LLU. (b) y]*" ~ 25, LLP. 

Figure 16: Detail of grid distribution, RAE 2822, x/c ~ 0.25 


Wall Integration, y* « 1 

Wall function, y* ~ 10, LLU 



Wall Integration, y* - 1 

Wall function, y* « 10, LLU 



(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 17: Solution residual, continuity equation, RAE 2822, Case 9 
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Wall integration, y* ^ 1 . 

Wall function, y* = 10., LLU 
Wall function, y* = 20., LLU 
Wall function, y* »= 40., LLU 
Wall integration (NPARC), y* = 2.5 



Wall Integration, y* 1 . 

Wall function, y* »= 10., LLU 
Wall function, y* »= 20., LLU 
Wall function, y* »= 40., LLU 
Wall integration (NPARC), y* = 2.5 



(a) SA turbulence model. (b) SST turbulence model. 

Figure 18: Upper and lower wall spacing distributions, RAE 2822, Case 9 


o Data, Cook et al. 

Wall integration, y^ « 1 . 

Wall function, y* »= 1 0., LLU 

Wall function, y* »= 20., LLU 

Wall function, y* = 40., LLU 

— Wall integration (NPARC), yl = 2.5 



o Data, Cook et al. 

Wall integration, y^ » 1 . 

Wall function, y* »= 1 0., LLU 

Wall function, y* »= 20., LLU 

Wall function, y* »= 40., LLU 

— Wall integration (NPARC), yl » 2.5 



X / C 


X / c 


(a) SA turbulence model. 


(b) SST turbulence model. 


Figure 19: Pressure coefficient distributions, RAE 2822, Case 9 
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Data, Cook et al. 

Wall integration, y* » 1 . 
Wall function, y* = 10., LLU 



Data, Cook et al. 

Wall integration, y* 1 . 
Wall function, y* »= 10., LLU 



(a) SA turbulence model. (b) SST turbulence model. 

Figure 20: Skin friction coefficient distributions, RAE 2822, Case 9 


Table 4: Computed forces and moments, RAE 2822, SA turbulence model. 


Grid 


total 

Lift 

pressure 

viscous 

total 

Drag 

pressure 

viscous 

Pitching 

moment 

Wall integration 

1. 

0.791 

0.791 

0.000 

0.0173 

0.0117 

0.0056 

-0.090 

Wall function (LLU) 

10. 

0.781 

0.781 

0.000 

0.0171 

0.0114 

0.0057 

-0.088 

Wall function (LLU) 

20. 

0.782 

0.782 

0.000 

0.0171 

0.0115 

0.0056 

-0.089 

Wall function (LLU) 

40. 

0.781 

0.781 

0.000 

0.0172 

0.0116 

0.0056 

-0.089 

Wall integration (NPARC) 

2.5 

0.784 

0.784 

0.000 

0.0175 

0.0118 

0.0057 

-0.089 

Wall function (LLP) 

25. 

0.769 

0.769 

0.000 

0.0169 

0.0114 

0.0055 

-0.086 


Table 5: Computed forces and moments, RAE 2822, SST turbulence model. 


Grid 


total 

Lift 

pressure 

viscous 

total 

Drag 

pressure 

viscous 

Pitching 

moment 

Wall integration 

1. 

0.768 

0.768 

0.000 

0.0164 

0.0109 

0.0055 

-0.086 

Wall function (LLU) 

10. 

0.764 

0.764 

0.000 

0.0162 

0.0108 

0.0054 

-0.085 

Wall function (LLU) 

20. 

0.756 

0.756 

0.000 

0.0161 

0.0107 

0.0054 

-0.084 

Wall function, (LLU) 

40. 

0.752 

0.752 

0.000 

0.0162 

0.0107 

0.0055 

-0.083 

Wall integration (NPARC) 

2.5 

0.763 

0.763 

0.000 

0.0163 

0.0110 

0.0053 

-0.084 

Wall function (LLP) 

25. 

0.742 

0.742 

0.000 

0.0159 

0.0106 

0.0053 

-0.080 
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D. ONERA M6 wing 

The fourth configuration is the ONERA M6 wing.^^ This geometry is well known and often used in CFD 
validation studies. The M6 wing is challenging not only for the highly three-dimensional nature of the 
flow, but the upper surface shock and degree of flow separation. These features are quite prominent even at 
3 and 4 degrees angle of attack. 

The unstructured grids for the M6 were generated with the VGRID v4.0 advancing-layer and advancing- 
front grid generation software package.^® VGRID has two types of spacing requirements: the “inviscid” 
spacing distributions are used in the advancing-front region of the grid, and the “viscous” spacing distribu- 
tions are used near the solid body in the advancing-layer regions of the grid where high stretching is required. 
The clustering of points normal to the surface was computed according to the VGRID stretching function,®® 
Sn = <5i[I -I- ri(I -I- r 2 )"'“^]"“^, where is the normal spacing of the layer, <5i is the spacing of the first 
layer, and the factors iq and iq are constants that determine the rate of stretching. (Note if iq is zero the 
stretching is geometric.) 


Table 6: Grid parameters for the ONERA M6 wing. 


Level 

Nodes 

ri 

f2 


Wall integration 

6.6x10® 

0.15 

0.02 

- 1 . 

Wall function, LLU 

6.1x10® 

0.10 

0.02 

~25. 

Wall function, LLU 

5.7x10® 

0.10 

0.02 

-50. 


Three different grids were generated for this study: a wall integration grid with a nominal y)*" of 1, and 
two LLU wall function grids with nominal yj*" of 25 and 50, and are listed in Table All of the grids use 
the same set of inviscid spacing distributions, and therefore, the surface grid for each grid has essentially 
the same surface spacing distributions. For the integration to the wall grid, the grid normal spacing is set 
such that an average normalized coordinate y®" is less than one for the first grid cell at the wall for the 
majority of the wing. The advancing-layer stretching rates, ri and r 2 , are 0.15 and 0.02, respectively. The 
total number of grid nodes is 6.6 million. For the second grid, the grid normal spacing is increased to 25, 
and the advancing-layer stretching rates, ri and r 2 , are 0.10 and 0.02, respectively. The total number of grid 
nodes for the second grid is 6.1 million. For the third grid, the grid normal spacing is increased to 50, and 
the advancing-layer stretching rates, iq and iq, are 0.10 and 0.02, respectively. The total number of grid 
nodes is 5.7 million. 

The grids for the isolated wing extend approximately 33 MAG away from the surface to a square outer 
boundary with a symmetry condition on the centerline. The maximum spacing at the wing grid outer bound- 
ary is approximately 0.05 MAG. The characteristic chordwise spacing at the leading edge is approximately 
0.001 MAG. The characteristic chordwise spacing at the trailing edge is approximately 0.002 MAG. There 
is no spanwise grid stretching at the leading or trailing edge. The characteristic spanwise spacings at the 
root and tip are approximately 0.01 and 0.004 MAG, respectively. A cylindrical volume of refinement is set 
around the wing surface such that within this volume the maximum grid spacing cannot exceed 0.02 MAG 
near the root and 0.005 MAG at the tip. This cylindrical volume of refinement was achieved by utilizing 
VGRID volume sources.®® The total number of wing surface points is approximately 107,000 for each of the 
grids. 


A comparison of the relative density of the wall integration grid and the LLU yj^ ^ 50 grid is shown in 
Fig. 1^ As expected, the near-wall grid density is significantly reduced with the larger wall spacing. 

The M6 wing was run at a Reynolds number of 14.0 million/chord (21.7 million/meter), freestream Mach 
number of 0.84, a freestream temperature of 540. OK, and at 3.06 and 4.08 degrees angle of attack. Solution 
residuals for the continuity equation are plotted in Fig. |22[ Only six orders of magnitude reduction in 
continuity equation residual for the SA calculations at a = 3.06° was achieved. Eight orders of magnitude 
reduction, or better, in the continuity equation was achieved with the SST model calculations, as well as, 
the a = 4.08° calculations. The Venkat limiter^® was required for the a = 4.08° calculations. The limiter 
was frozen after the first 4000 iterations, which improved the residual convergence after that point. 

Figure 23 is representative of the difference between the various simulations and experimental data. 
Very good agreement of pressure coefficient distributions was obtained between wall integration and the two 
wall function grid solutions at a = 3.06° degrees, though neither the wall integration nor the wall function 
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Figure 21: Detail of slice through the grid, x/c ^ 0.374, ry = 0.25, ONERA M6 wing. 


calculations accurately placed the location of the shock on the wing. Relatively good agreement of pressure 
coefhcient distributions occurred at a = 4.08° degrees, though the shock location moved slightly upstream 
with the wall function grids. 

Surface streamlines for the wall integration, y]*" « 1, and wall function, y^ « 25, at a = 4.08° degrees are 
plotted in Fig. Strong cross-flow occurs on the outboard region of the wing, as seen by the streamlines 
moving outboard, toward the wing tip, downstream of the upper surface shock. The SA wall model on 
the y^ ~ 25 grid, to some degree, mimicked the wall integration streamlines. The SST wall model on 
the same grid did not show nearly the outward turning. The surface streamlines are calculated from the 
average velocities in the elements in the grid attached to the wing boundary. Wall integration grids, with 
wall spacing of 5 wall units or less, are assumed to provide a reasonable representation of the surface flow. 
It is questionable whether wall function grids can, in general, provide the correct average velocity held. 
The further the boundary element extents beyond the viscous sub-layer, the less linear velocity profile is. 
Therefore the surface streamlines from wall function grids have some inherent error. 

A summary of the lift, drag and pitching moments are listed in Tables[^and[^ The pressure contribution 
to the lift and drag were fairly consistent between solutions obtained on the wall integration and the wall 
functions for both turbulence models and angles of attack. The viscous contribution to the drag computed 
by the wall function on the y^ ~ 50 grid, though, was 40% below that of the wall integration viscous drag 
using the SA wall function model for both angles of attack. The SST wall function model computed similar 
viscous drag numbers as the SA model. 


V. Summary 

The implementation of two wall function models into the node-based unstructured-grid CFD code Fun3D 
has been described. The hybrid wall models of Knopp et ah,® were used to simulate flow along a zero pressure 
gradient flat plate, an asymmetric diffuser, a 2D airfoil and a 3D semi-span wing. Reasonable agreement 
with canonical and experimental data, comparing pressure and skin friction, was observed for the flat plate 
and asymmetric diffuser test cases. Computed friction velocity using log-layer packed type grids were more 
constant across a range of initial wall spacing for the flat plate. Log-layer unpacked type grids produced 
velocity, skin friction and pressure coefficient distributions very similar to those obtained with log-layer 
packed or wall integration grids for the RAE airfoil. The pressure distributions were close when comparing 
calculations using wall function grids and using the wall integration grids. The shock location and spanwise 
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Wall Integration, « 1 . Wall Integration, y* = 1 




(a) SA turbulence model, a = 3.06°. 


(b) SST turbulence model, a = 3.06°. 


Wall Integration, y* 1 . 
Wall function, LLU, y!j^ = 25. 



(c) SA turbulence model, a = 4.08°. 


Wall Integration, y* 1 . 
Wall function, LLU, y* = 25. 



Figure 22: Solution residual, continuity equation, ONERA M6 wing. 


Table 7: Computed forces and moments, ONERA M6, SA turbulence model. 


Grid 


a 

deg. 

total 

Lift 

pressure 

viscous 

total 

Drag 

pressure 

viscous 

Pitching 

moment 

Wall integration 

1. 

3.06 

0.260 

0.260 

0.000 

0.0165 

0.0116 

0.0049 

-0.180 

Wall function (LLU) 

25. 

3.06 

0.263 

0.263 

0.000 

0.0155 

0.0118 

0.0037 

-0.184 

Wall function (LLU) 

50. 

3.06 

0.264 

0.264 

0.000 

0.0149 

0.0118 

0.0031 

-0.184 

Wall integration 

1. 

4.08 

0.359 

0.359 

0.000 

0.0266 

0.0219 

0.0047 

-0.256 

Wall function (LLU) 

25. 

4.08 

0.354 

0.354 

0.000 

0.0252 

0.0218 

0.0034 

-0.251 

Wall function (LLU) 

50. 

4.08 

0.354 

0.354 

0.000 

0.0246 

0.0219 

0.0027 

-0.251 
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Data, Schmitt et al. 



Data, Schmitt et al. 
Wall Integration, yif = 1 . 



25 . 

50 . 


x/c 


x/c 


(a) SA turbulence model, a = 3.06°. 


(b) SST turbulence model, a = 3.06°. 


o Data, Schmitt et al. 


Wall Integration, y* ^ 1 . 



o Data, Schmitt et al. 


Wall integration, y:^ « 1 . 



x/c 


x/c 


(c) SA turbulence model, a = 4.08°. 


(d) SST turbulence model, a = 4.08°. 


Figure 23: ONERA M6 wing, j] = 0.80. 


Table 8: Computed forces and moments, ONERA M6, SST turbulence model. 


Grid 


a 

deg. 

total 

Lift 

pressure 

viscous 

total 

Drag 

pressure 

viscous 

Pitching 

moment 

Wall integration 

1. 

3.06 

0.264 

0.264 

0.000 

0.0166 

0.0116 

0.0050 

-0.184 

Wall function (LLU) 

25. 

3.06 

0.264 

0.264 

0.000 

0.0153 

0.0117 

0.0034 

-0.185 

Wall function (LLU) 

50. 

3.06 

0.264 

0.264 

0.000 

0.0147 

0.0117 

0.0030 

-0.184 

Wall integration 

1. 

4.08 

0.364 

0.364 

0.000 

0.0268 

0.0220 

0.0048 

-0.261 

Wall function (LLU) 

25. 

4.08 

0.353 

0.353 

0.000 

0.0252 

0.0218 

0.0034 

-0.249 

Wall function (LLU) 

50. 

4.08 

0.353 

0.353 

0.000 

0.0247 

0.0218 

0.0029 

-0.250 
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(c) SA turbulence model, ^ 25. (d) SST turbulence model, 25. 

Figure 24: Upper surface streamlines, ONERA M6 wing, a = 4.08°. 
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flow features on the ONERA M6 wing were more difficult to calculate with the wall function models. The 
integrated viscous forces were off by 40% for the « 50 grid. 
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